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Abstract 

In contemporary image and vision analysis, stochastic approaches demonstrate great flexi- 
bility in representing and modeling complex phenomena, while variational-PDE methods gain 
enormous computational advantages over Monte-Carlo or other stochastic algorithms. In com- 
bination, the two can lead to much more powerful novel models and efficient algorithms. In 
the current work, we propose a stochastic-variational model for soft (or fuzzy) Mumford-Shah 
segmentation of mixture image patterns. Unlike the classical hard Mumford-Shah segmentation, 
the new model allows each pixel to belong to each image pattern with some probability. We 
show that soft segmentation leads to hard segmentation, and hence is more general. The mod- 
eling procedure, mathematical analysis, and computational implementation of the new model 
are explored in detail, and numerical examples of synthetic and natural images are presented. 

Keywords: Segmentation, soft, mixture, pattern, ownership, probability simplex, Modica-Mortola, phase- 
field, symmetry, supervision, Egorov's theorem, Poincare inequality, existence, AM algorithm. 
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1 Introduction: Soft vs. Hard Segmentation 



Segmentation is the key step towards high-level vision modeling and analysis, including object 
characterization, detection, and classification. There have been some recent developments indicat- 
ing that certain high-level visual tasks such as global scene interpretation might be able to bypass 
segmentation |24( I23j . Nevertheless, segmentation still remains perhaps the most important and 
inspiring task to date in low- or middle-level vision analysis and image processing. 

The segmentation problem can be formulated as follows. Given an image / on a 2-dimensional 
(2D) domain Q (assumed to be bounded, smooth, and open), one seeks out a closed "edge set" F, 
and all the connected components Oi, . . . , ^Ik of O \ F, such that by certain suitable visual measure 
(e.g., textural or photometric), the image / is discontinuous along F while smooth or homogeneous 
on each segment Qi. Each image patch Jj = is also called a pattern, and J7j its support. 

We shall call this most common practice "hard" segmentation. A hard segmentation partitions 
the image domain Q along a definitive edge set F, and outputs non- overlapping pattern supports 

r^i, . . . , 

The present work introduces the notion of "soft" segmentation. Mathematically, a hard seg- 
mentation amounts to the partition of the unit using indicator functions: 

K 

1^(2:) = ^ ln,(x), a.e. x = {xi,X2)eQ. 
1=1 

A soft segmentation seeks out instead a softer partition of the unit: 

K 

ln{x) = Y,P^i^)^ (1) 

i=l 

where pi's are continuous or smoother functions. Formally, each pi could be considered as the 
mollified version of 1q.{x). 

In the stochastic literature of image analysis and modeling, the above notion of soft segmentation 
is closely connected to mixture image models (e.g., ,2Qj). Suppose a given image / is composed from 
K unknown patterns: 

ca = 1, UJ = 2, . . . , UJ = K, 

where uj denotes the pattern label variable. At each pixel x G fi, uj{x) £ {I,-- - ,K} becomes a 
random variable. Then the pj's in carry the natural stochastic interpretation: 

Pi{x) = Prob(u;(x) = i), i = 1 : K. 

For this reason, each pi shall be called the ownership of pattern i. Instead of the repulsive ownership 
in a hard segmentation, a soft one allows each pattern to "own" a pixel with some likelihood. 

Soft segmentation is more general since it can lead to natural hard segmentation under the 
maximum likelihood (ML) principle. Given a soft segmentation {pi{x) : i = 1 : K}, one can define 
for each pixel x G J7 its unique owner u;*(x) by: 

uj^{x) = argmax^gi.^ pU^), (2) 

and if the maxima are non-unique, accept the largest index from the argmax pool. The segments 
are then defined by 

fli = uj^^{i) = {x G I uj^{x) = i}, i = 1 : K, (3) 
which leads to a natural hard segmentation. ((2)) and Q are called the hardening formulae. 
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Soft segmentation has been motivated by practical analysis of natural images. Patterns in 
natural scenes often do not have clear-cut boundaries. In Figure ^ for example, there does not 
seem to exist a "hard" boundary between the grass and sand areas. If one draws an oriented line 
as shown in the figure, it makes more sense to state that along the arrow the pattern transits 
from being "more" sand like to being "more" grass like. Such consideration favors the following 
stochastic view that along the arrow, the ownership 

Prob(a;(x) = grass) increases, while Prob(c<j(3;) = sand) decreases. 



Figure too large for arXiv 
For the figure and the complete manuscript, 

go to: www.math.umn.edu/~jhshen 



Figure 1: Natural images often do not have clear-cut "hard" boundaries between different patterns. 
Along the arrow, for example, one only observes that the sand pattern gradually becomes a grass 
pattern. Such a "soft" view is the stochastic view on the segmentation problem. 

In the present work, we propose a new stochastic-variational soft-segmentation model for the 
following celebrated Mumford-Shah model |3H I33j : 

mmE[u,T\I]=n\T) + a [ \Vu\^ + \ f {u - if , (4) 

where stands for the ID Hausdorff measure, which is simply the length when F is regular 
enough. For notational conciseness, the default area-element symbol dx = dxidx2 will he omitted 
in most integral formulae. 

As stated in the abstract, the stochastic softness induces more flexibility and universality in 
modeling, while the variational-PDE approach facilitates rigorous mathematical analysis as well as 
more efficient computational implementations compared with purely stochastic approaches includ- 
ing, e.g., the Monte-Carlo method or Gibbs' sampling |ll [TCI US] . 

The paper has been organized as follows. Section 2 builds up the soft Mumford-Shah (SMS) 
model under the Bayesian rationale and the MAP estimator \l\V I32j , which are the formal stochas- 
tic foundations of the present model. In Section 3, the prior energy on the ownerships piS is 
developed based on the celebrated work of Modica and Mortola ^30] on phase-field modeling and 
F-convergence approximation in material sciences and phase transitions. In Section 4, we analyze 
the main mathematical properties of the proposed SMS model, including the admissible space, 
hidden symmetry and symmetry breaking via weak supervision, and the existence theorems. In 
Section 5, we then derive the system of Euler-Lagrange equations of the SMS model for which 
the role of the probability simplex constraint is discussed in detail. Section 5 also introduces the 
alternating-minimization algorithm to compute the Euler-Lagrange equations. Finally, the numer- 
ical performance of the SMS model is demonstrated in Section 6 via both synthetic and natural 
test images that are sufficiently representative and generic. 
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Throughout the manuscript, the notation E\X, y | Z] in the deterministic setting always de- 
notes a quantity (often a functional) E that depends on X, Y, and Z but with Z given or fixed. 
Similarly, E[X \ Y^Z] still denotes E[X,Y \ Z] modulo some additive quantity g[Y, Z] that is of- 
ten unimportant as far as the optimization on X (given Y and Z) is concerned. These notations 
therefore have been closely inspired by conditional probabilities in the stochastic setting. 

2 Bayesian Rationale to the New Model and Gaussian Mixture 
2.1 Bayesian Rationale 

Segmentation can be done in some feature spaces such as gradient-like highpass features or Gabor 
features (e.g., [A^ I44| I45j ) . The Mumford-Shah model easily extends to such general features 
(e.g., even though it was originally formulated only for intensity fields. For maximal clarity 

in exposing the core ideas of the current work, we shall also focus only on the latter, while leaving 
as canonical exercises to adapt the new model for any given feature distribution. 

Let K be the total number of intended patterns. As in [431 145j . K could also be treated as 
an unknown to be optimally estimated, which however does not add much to the most significant 
contribution (i.e., the modeling and computation of the "softening" procedure) of the present work. 

Given an image input / = I{x) on a bounded, regular, and open domain fi, the primary goal 
of soft segmentation is to compute the ownerships: 

Pi{x),P2{x),...,pk{x). 

Define P{x) = {pi{x),p2{x), . . . ,pKix)), and 

= convex hull of ei, . . . , eV, 

where the (cj | i = 1 : K) denotes the canonical Cartesian basis of . A/^_i is often called the 
canonical [K — l)-simplex, or the probability-simi)lex in M^. Then 

P -.n^ Ak-1, x^P{x), 

meaning that the total ownerships always add up to 100% at any pixel x G 0. 

Associated with each pattern label a; = i is a smooth function Ui{x) £ H^(Q), similar to the 
original Mumford-Shah model. Here the Sobolev space H^{Q) is defined by [SU] 

H^{n) = {u£ L^{n) I Vn e L2(S7,IR2)}_ 

Define U{x) = {ui{x),U2{x), . . . ,UKix)). Then the goal of soft segmentation is to estimate the 
optimal vectorial pair of ownerships and patterns given an image /: 

(P*, U*) = argmax(p [7)Prob(P, U \ I). 

By the Bayesian formula j2H l,S2j. the posterior given / is expressible via 

Prob(P,(7 I /) = Prob(/ | P, J7)Prob(P)Prob(?7)/Prob(/), 

assuming that the mixture patterns U and the mixture rules P are independent. We shall call the 
first term a "mixture generation" model, since it reveals how the image data should look like given 
the information of the patterns and their ownerships. 
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By taking the logarithmic likelihood E[-] = — logProb(-), or formally Gibbs' energy in statistical 
mechanics [^1^], one attains the soft segmentation model in its "energy" form: 

min^ E[P, U\I]=E[I\ P, U] + E[P] + E[U], (5) 

modulo an insignificant additive constant. 

Assuming that all the pattern channels are independent from each other, one has 

K 

E[U] = E[ui,. . . ,uk] =Y,^i''i I 

1=1 

For Sobolev-regular patterns, one may impose the homogeneous Sobolev energies: 

E[ui I i] = E[ui] =a iVujp, i = l:K, (6) 
Jn 

for some scalar weight a that models the visual sensitivity to intensity roughness. Unlike the original 
Mumford-Shah model, the energy for each channel has been defined on the entire image domain 
Q instead of on each "hard-cut" patch fij. Thus the energy form @ must carry out extrapolation 
for practical applications. Long-range extrapolations are, however, often unimportant after being 
weighed down by their negligible ownerships pi^s. 



2.2 Gaussian Mixture with Smooth Mean Fields 

In this subsection we discuss the mixture generation model Prob(/ | P, U) or E[I \ P,U]. 

Assume that the patterns are all Gaussian's with mean fields ui,U2, . . . , uk- For simplicity also 
assume that they share the same variance cr^ (which readily generalizes to the more general case 
with variations). Then at any given pixel x £ Q, 

{I I u;{x) = i)r^ N{ui{x), 0-2), i = l:K. 

Define the Gaussian probability density function (p.d.f) 

(T\ \ 1 ( {I-^y 

g[i \ m,a) = exp 



V2^a V 

The the p.d.f of the mixture image / at any pixel x is given by 

K K 

Prob(/(x) I P{x),U{x)) = ^Prob(/ | cj(x) = i) Prob(t(;(x) = i) = ''^g{I \ Ui{x),a)pi{x). 

1=1 1=1 

Thus ideally the "energy" for the mixture generation model should be given by 

E[I \ P,U] = — /u / log y^S^^ I '^ii^) ^ ^)Pii^) ' fo'^ some scalar ^ > 0, 

\i=i J 

provided that, given two fields P and U on i}, for any two distinct pixels x and y, 

{I{x) I P, U) is independent of {I{y) \ P,U). 



(7) 



5 



In the current work, we shall adopt a reduced form of the complex formula , which is simpler 
and easier to manage both in theory and for computation. Assume that each soft ownership Pi{x) 
is closer to a hard one Pi{x) ~ l^^(x) for i = 1 : K. Then 



K \ / K 



log [^2^^^ ^ Ui{x),a)pi{x) j ~ - log (^^S-U I Ui{x),a)ln^{x)j 

K 

'^log g{I \ui{x),a)l^^{x) (a.e.) 

i=l 
K 

^ log 5(1 I Ui{x),a)pi{x) 



^i=l / \i=l 

K 



1=1 

K 



^ ~ ^ij(a;))^K(a^) + const., 



2a2 . 

1=1 

where the additive constant only depends on cr and K. This suggests the following convenient 
energy form for the mixture generation model: 



E[I \P,U] = \j^ [y1{I - Ui{x)fpi{x)^ , 



(8) 



which amounts to a weighted least-square energy j42j . The weight A reflects visual sensitivity to 
synthesis errors. 

In combination of (jSJ, and (jSJ, the new soft segmentation model takes the form of mini- 
mizing 

K K 

E[P,U\I] = \y^ / [I -u,{x)fp,{x) + ay^ / |V7/,|2 + S[P]. (9) 

i=l ^=l 

Notice that here the ownership distribution P "softens" the "hard" segmentation boundary F in 
the original Mumford-Shah model To complete the modeling process, it suffices to properly 
define the prior or regularity energy E[P\, which is the main task of the next section. 

3 Modica-Mortola's Phase-Field Model for Ownership Energy 

To generalize but not to deviate too far from classical hard segmentation, it is natural to impose 
the following two constraints: 

(a) each pattern ownership Pi{x) has almost only two phases: on (corresponding to Pi = 1) and off 
(to Pi = 0), and the transition band in between is narrow; and 

(b) the soft boundaries, or equivalently the transition bands, are regular, instead of being zigzag. 

In combination, one imposes the following Modica-Mortola type of energy with a double-well po- 
tential ISni: Pi G H^{^), 

i?,N= j[^(9e|Vp.p + Mi^^), i = l:K. (10) 

Here e ^ 1 controls the transition bandwidth. Since e <C 1, the second term necessarily demands 
Pi ~ or 1 to lower the energy, which well resonates with the expectation in (a). The first term, 
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weighted by the small parameter e, amounts to a regularity condition on each pi, which meets the 
requirement in (b). 

Energies in the form of (jlUj) are very common in material sciences, including the theories of 
liquid crystals and phase transitions jT^ EI ■ Mathematically they have been well studied in the 
framework of F-convergence |20i, which we now give a brief introduction in the present context. 
We also refer the reader to the works of Ambrosio and Tortorelli |2j on the F-convergence 
approximation to the classical Mumford-Shah segmentation model. 

Recall that for any q{x) S L^{Q), its total variation as a Radon measure is defined by [7| I18j 



TY[q] = [ \Dq\ = sup (g, V • g 



g&C},{fl,B2) 

where stands for the unit disk centered at the origin in M?. (The TV measure was first introduced 
into image processing by Rudin, Osher, and Fatemi j34j.^ Define for any q E L^(J7), 



Eo[q] 



{TV[g], if g = or 1, a.e. on Q, 

oo, otherwise. 



As a result, a finite energy Eolq] necessarily implies that q has two phases only, and E^lq] = TY[q] = 
Peic{q~^{l)) is the perimeter of the support region V = q~^{l). 
Further define 

Lfo_^](17) = {qe L\Q) I q{x) G [0, 1], V a; G 

to be a subspace of L^(i7) (as a metric space). Then Modica and Mortola's well known results 
in readily leads to the following theorem. 

Theorem 1 (F- Convergence Approximation of a Two-Phase TV) For any q G L^q-^^{Q) \ 
H^{Q), we extend the definition of Eel-] GU) defining Eglq] = +oo. Then 

Egr Eq in the sense of T- convergence in the metric space LjQ-^j(i7). 

That is 

(i) for any q^ ^ q in Ljq -^j(f7) as e ^ 0, 

lim inf E^ [g^] > Eq [q\ ; and 

(ii) for any q G Ljq there exists some sequence {q* \ e), such that q* q as e ^ 0, and 

Ihn Ee[q;]=Eo[q]. 

e— >0 

We refer the reader to Modica and Mortola [HO] for a proof (with some necessary modification). 
Here we only point out that the "tight" sequence {q* \ e) in (ii) can be constructed using a smooth 
sigmoid transition across the hard boundary of a given two-phase function q. Recall as in the theory 
of neural networks [S] that a sigmoid transition between and 1 is achieved by 

1 

a(t) = 7, — oo < t < oo. 

1 -I- e * 

The scaling parameter e participates in the transition by the form of a{t/{3e)). In particular, e 
indeed corresponds to the width of the transition band when t is a distance function. 

This theorem reveals the close connection of the particular choice of E^ \pi] in (|1L)|1 to the original 
Mumford-Shah model. 
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Proposition 1 Suppose p^'s "optimally" (i.e., by the above sigmoidal transition) converge to a 
given 2-phase pattern ly(a;) with a regular hard boundary T = dV . Then, 

E,\p,] ^ length(r) = / \Dlv{x)\ . 

Jn 

Similar results have appeared in the earlier influential works of Ambrosio and Tortorelli |^ [2] 
on the F-convergence approximation to the Mumford-Shah model. The technique has also been 
extensively applied in image computation and modeling ^1 UHl OH OHl BHl to overcome the 
difficulty in representing and computing the free boundary T. 

To summarize this section, we propose the following energy model for the ownership distribution 
P{x) = {pi{x),p2{x), . . . ,pk{x)): 

E,[p] =f:E,[p.] = jz [ + 

.=1 i=i -^^ V e 

One, however, must realize that different ownerships are not decoupled by this energy though it has 
appeared so. The energy -Ee[-P] must be coupled with the constraint of the probability-simplex: 

K 

P : Q ^ Ak~i, or = 1, Pi >0, \/ x £ Q. 

i=l 

In particular, for small e, although implies that each ownership pi tends to polarize to or 
1 independently, they have to cooperate with each other under the above simplex constraint to 
optimally share the ownerships. 




4 Soft Mumford-Shah (SMS) Segmentation 
4.1 The Model and Admission Space 

Combining the preceding two sections, we have developed the complete formula for soft Mumford- 
Shah segmentation with iT-patterns: 

mmE[P,U\I] = Xy2 / - + « V / iVn^l^ + V / f 9e| Vp.|^ + ^^'^^ ~ ^^^^^ V (12) 

with the constraint that 

P : 0, ^ Ak-1, the probability (K — l)-simplex, 

i.e., Pi > 0,i = 1 : K , and Yli^iPi — ^- discussed previously, it is this simplex constraint that 
induces coupling among different channels into the seemingly decoupled model (|12() . 

Besides the simplex constraint, the last term in the energy (|12() requires G H^{^) for i = 1 : K. 
Similarly, the second term requires each pattern m S H^{Q,). Then with the assumption that 

"the given image / G L^(J7)," 

E[P,U I /] is well defined and finite for any admissible patterns U and pattern ownership distri- 
bution P: 

admK = {{P,U)\pi,UieH^{n),i = l:K; P:Q^Ak-i}. (13) 
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4.2 Breaking the Hidden Symmetry via Weak Supervision 

Let Sk denote the permutation group of {!,... Each permutation a G Sk is a 1-to-l map: 

a:{l,...,K]^{l,...,K], 

so that (17(1), . . . , (y{Ky) is a re-arrangement of {1, ... , K}. For any i^-tuple F = (/i, . . . , fx), one 
defines 

Fu = (/o-(l) : /cr(2) ) • • • J fa-{K))- 

Theorem 2 (Hidden Symmetry of SMS) For any a G Sk, 

E[P^,U^\I]=E[P,U\I]. 

In particular, suppose 

{P*,U*) = argmin(p,t7)eadm^,ii;[P,[/ | /] 
is an optimal pair. Then for any a G Sk, (.P^, U*) is a niinimizer as well. 

The proof is straightforward and thus omitted. Such symmetry not only worsens the non- 
uniqueness of the non-convex energy functional in H12|) . but also potentially jitters intermediate 
solutions in iterative computational schemes (i.e., hysterical transitions in the admissible space). 

To break the permutation symmetry, we turn to a weak supervision scheme in which a user 
specifies K distinct domain patches: 

Qi,Q2, ■ ■ ■ , Qk, 
and imposes the symmetry- breaking conditions: 

P^Q, = '^ii' ^'-^ = 1: K, (14) 

where 6ij denotes Kronecker's delta. That is, a user requires each given patch Qi to be a "pure" 
pattern exclusively labelled by i. Computationally this weak supervision process can be automated 
based on multiscale patch statistics as in the contemporary works on scene recognition |23( I24j , or 
more generally, the learning theory |1U1 141j . 



Figure too large for arXiv 
For the figure and the complete manuscript, 

go to: www.math.umn.edu/~jhshen 



Figure 2: Examples of a 3-phase supervision and a 4-phase supervision to break the symmetry in the 
model. Such weak supervision can also be automated based on multiscale patch statistics j24| iS^I . 
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4.3 Existence Theorems for Non-Supervision and Supervision 

In this subsection, we establish the existence theorems for the soft Mumford-Shah segmentation 
model H12|) without or with the supervision H14|) . 

Theorem 3 (Existence Theorem for Unsupervised SMS) Suppose I € L'^{Q). Then for 
any positive modeling parameters (A,Q,e), a minimizer to the unsupervised soft Mumford-Shah 
model HE) must exist. 



We will need the following lemma for the proof. 

Lemma 1 Let {fn \n) he a sequence of functions in L'^{Q), and (p" | n) a sequence of nonnegative 
measurable functions on and valued in [0, 1] . Suppose 

(i) p* , a.e. on and JqP* > 0; and 

(ii) fnP^ — ^ /^'^ some A> and n = 1 : oo. 
Then there exists some function p G L'^{Q), such that 

(a) p > and f^p = 1, and 

(b) for some fixed B > 0, < B for n = 1 : oo. 

Proof. Denote the Lebesgue measure of a measurable set W by \ W\. Since p* > and JqP* > 0, 
there must exist some c > 0, such that 

V = {x €z n \ p* > 2c} has a finite but positive measure. 

On the other hand, by Egorov's theorem ^1] on a.e. convergence, there must exist a subset W C V, 
such that 

(a') \V-W\< ^, and hence \W\ > 0, and 
(b') p* uniformly on W. 

In particular, there exists some N, such that for any n > N, p"^ > c on W. Define 



p{x) = e L^m. 



Then J^p = 1, and for any n > N, 



n 



c\W\ Jw - c\W\ Jn'''^ - c\W\ 
Thus by the Schwarz inequality (or E[X]^ < E[X^] in probability theory), 

/..pls(X/^)'^(^)'^ n>N 

The lemma holds if one defines B = max ^(^^j]^) ) \Jq fip\ > • • • > |/n fNp\ 1 • D 
We are ready to prove Theorem |21 
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Proof. Take the special pattern distribution: 

Ui = 0, i = l : K; pi = 1, Pj =0, j = 2 : K. 

Then 



E[P, U\I] = \ [ /2 < oo. 
Jn 



Thus the infimum of the energy must be finite. Let (-P", [/" | n) C admi^ (see (|13j) ) be a minimizing 
sequence for the soft Mumford-Shah energy (|12j) . 

Due to the third term in the energy and the simplex constraint, for each channel i, (p" | n) must 
be bounded in H^{Q). By the L'^-weak compactness, there must exist some P* S L^(J7,M^), and a 
subsequence of (-P" | n), which after relabelling shall still be denoted by (-P" | n) for convenience, 
such that 

pn _^ p* in M^), n ^ oo. 

Then by the lower semi-continuity of Sobolev measures, 

9e [ |Vp*p < liminf 9e [ jVp^p, i = l:K. (15) 
Jn Jn 

Furthermore, with possibly another round of subsequence refinement, one can assume 

P"(x)^P*(x), a.e. xGQ, n^oo. 
Since the probability simplex Afc_i is closed and P"(x) G Ak-i, one concludes that 

P*{x) £ Ak-1, a.e. x £ Q. 
And by Fatou's Lemma jl4l I25j . one has 

'"^'^-^-'''Sliminf /•'"?" -^?»\ .^l:K. (16) 



(In fact, the equality holds by Lebesgue's Dominated Convergence pSj.l 

After the above subsequence selection on /'"'s, one naturally has an associated subsequence 
of ([/" I n), which for convenience is still denoted by (U"^ \ n) after relabelling. For each specific 
channel i, we then consider two scenarios separately. 

Suppose p*{x) = 0, a.e. x £ Q. We then define for that channel 

u*{x) = 0, x£n. (17) 

Such a channel is called a "dumb" channel. 

Otherwise, one must have JqP* > 0, and from the first term in H12|l . 

/ (ti" — I)'^Pi < const., n = 1 : oo. 
Jq 

Since J^I^pf < /n-^^' triangle inequality, 

/ {'ai)'^Pi < const., n = 1 : cx), 
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where the constant only depends on / and the model parameters. Then by Lemma ^ there exists 
some Pi{x) > 0, with J^Pi = 1, some constant Bi > such that 



UiPi < Bi, n = 1 : oo. 

On the other hand, by the second term in the energy ()12() . 

/ |V<|2 < Q = Q(/,A,a,e), n=l:cx), 
Jn 

for some constant Cj independent of n. Then by the generalized Pointcare inequality jl3ll25j on 

\\w - {w,Pi)\\L2 < Ai\\Vw\\L2, 

where Ai = Ai{pi,Q,) is independent of w G H^{Q), one concludes that 

Wu'^Wl^ < Di = Di{Ai,Bi,Ci), n=l:(X), 

for some constant Dj. As a result, {uf \ n) must be bounded in By the L^-weak compactness 

of bounded ^-sequences, there is a subsequence of (ti" | n), for convenience still denoted by (n" | n) 
after relabelling, such that 

^ n* G n ^ oo, 

converging in the sense of both LP' and almost everywhere. Then by the lower semi-continuity, 

|Vu*P < liminf |Vm"P. 



(18) 



Finally, since u^{x) and pf{x) Pi{x), a.e. x £ Fatou's Lemma gives 

/ K-I)V<liminf / K-I)^p:. 
Combining both cases just analyzed above, we have established that 

K K K K 

i=l i=l n^oo 

Together with (|15() and (|16|) . this implies 

E[P*,U* I /] < liminf£;[P",J7" | /] = inf E[P,U \ I], 

n^oo {P,U) 

and hence (P*, U*) must be a minimizer. □ 
The proof reveals an important behavior of the model (|12l) . If certain channel i becomes dumb 
(i.e., pI = 0), it has often been introduced unnecessarily in the first place, and the associated 
optimal pattern u* could be any featureless constant image. 

Theorem 4 Suppose I G L^($7). Then an optimal pattern- ownership pair must exist to the soft 
Mumford-Shah segmentation model with supervision assuming that each patch Qi has a 

positive Lehesgue measure \Qi\ > 0. 

Proof. The proof is almost identical to the unsupervised case above, and simplifies substantially 
by noticing that no channel could become dumb due to supervision. Furthermore, the functions 
Pi's in the above proof can be directly set to be 

1 ^ ^ 

without the necessity of turning to Lemma ^ □ 
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4.4 Mixture of Homogeneous Gaussians 

When each pattern z is a homogeneous Gaussian N{mi,a) with a distinct mean value rrii, one has 

Ui{x) = rrii, X ^ i = 1 : K. 
Define m = (mi, . . . jTUk)- As a result, the soft Mumford-Shah model simplifies to 

mm E[P,m\ I] =Xy [{I- m,fp, + y f fgelV^p + MLzM!^ (19) 

Theorem 5 Suppose I G L^{Q,). Then a minimizer pair {P*,Tn*) to E[P,m \ I] exists for both 
the unsupervised and supervised cases. 

The proof can be derived readily from the previous general cases and is hence left out. When 
K = 2, a similar model was proposed earlier by Shen 37 under the symmetrization transform: 

1 - Z{x) 1 + Z{x) r -, n 

Plix) = — P2{x) = — z e [-1, 1]. 

The model p9() could be considered as the soft version of Chan and Vese's model |S] from 
the point of view of region-based active contours. Chan and Vese have demonstrated that such a 
piecewise constant Mumford-Shah model (or the CV model as popularly referred to in the present 
literature) is already powerful enough for a number of applications including medical imaging. 

5 Euler-Lagrange Equations and Computation on (i^' — 1)-Simplex 

5.1 Euler-Lagrange Equations on (K — 1)-Simplex 

To minimize the energy for the soft Mumford-Shah segmentation 

E[P,U\I]=xf2 [{u^-lfp^ + ay [ iV^d' + E/ f9.|Vp.p + Mjl^y (20) 

one resorts to its gradient-descent flow or Euler-Lagrange equations. In this section, we discuss 
these equations and their practical computational schemes. 

The first-order partial variation on U given P leads to, for i = \ : K, 

aAui + A(/ — Ui)pi = 0, on 17; -- — = 0, along dQ, 

on 

where n stands for the outer normal vector field along dO,. Thus the Euler-Lagrange equations on 
the patterns are all in the form of linear Poisson equations with variable coefficient fields: 

-aAui + {\pi)ui = fi, i = l : K, 

with Neumann adiabatic boundary conditions, where the source terms are fi{x) = Xpi{x)I{x). 

The first-order variation on the ownerships P is carried out on the probability {K — l)-simplex 
Ak-i, which is a compact manifold (with border) of codimension 1 embedded in M^. Chan and 
Shen [HI developed a general framework for modelling and computing image features that "live" on 
general manifolds, and especially those that are embedded in M^. We shall follow the approach 
there. 
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Without the simplex constraint on the ownerships, for any given U , the first order variation of 
the soft energy E under P — > P + 5P is given by 



K K 

6E = S2Vi6pidx+ / y^Vidpidn^, 
Jn Jdn 



where Tl is the 1-D Hausdorff measure along and 

Vi = X{ui - If - ISeApi + 2e-^pi{l - pi){l - 2p,), (21) 

Vi = l8s^, along aJl. (22) 
on 

Define V = {Vi, . . . , Vk) and v = {vi, . . . , vk)- Then 

6E= [ V-5Pdx+ [ v-6Pdn^, 
Jn Jdn 

which holds for any free variation of P in M^, or one writes in the free-gradient form 

vL + v 



In reality, P G /S.k-i- Let Tp/S.K-1 denote the tangent space of A/^_i at any single point 
P G A A', and 

TT : TpM.^'^ Tp Aa'-i 

the orthogonal projection onto the tangent space in M^. Since the normal direction of the tangent 
plane is given by \k/^/K = (1, . . . , l)/\/X, the projection operator is explicitly given by, for any 

1 ^ 

■k{w) = w — 1k{w^1k) / K = w — {w)1k, with {w) = — Wj- 



The constrained gradient of E on is therefore given by 

dE ( dE 



1=1 



In particular, the system of Euler-Lagrange equations on P given U is given by 

(Vi{x)-{v){x) = o, x£n, 

\vi{z) - {v){z) =0, z£dn, 
for i = 1 : K. The coupling among different channels is evident from these two formulae. 



(23) 



Lemma 2 Suppose P : 0, ^ A^'-i- Then for any z £ d^l, {v){z) = 0, where the boundary "flux" 
V is defined in H2I^) . 
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Proof. This is obtained by direct computation: at any z S dQ, 

1 18e dpi 

i=l i=l 

_18e_5_/^ \ _ 18£ _ 

□ 

As a result, the boundary conditions in ()23p simphfy to the ordinary Neumann conditions 
dpi/dn = Q,i = \ : K. Combining all the above derivations, we have established the following 
theorem. 

Theorem 6 (Euler-Lagrange Equations) The system of Euler- Lagrange equations ofE[P,U \ 
I] are given by 

-aAui + {Xpi)ui = {\pi)I, 

(24) 

-lSe/^Pi + 2e-^Pi{l-pi){l-2pi) = {V) -\{ui-lf, i = 1 : K, 

on 0,, all with Neumann boundary conditions along dO,. Here V = V{P,U) is defined as in \21)). 
Furthermore, under supervision J j^| ), the ownerships must satisfy the interpolation conditions: 

Or equivalently, the equations on pi 's in ^ hold onQ \ (uf^^Qi) with 



Neumann 



conditions along dQ, and Dirichlet conditions along Ll^^dQi : Pilgq. = S-, 



Similarly, one has the following result for the piecewise constant SMS model H19|) . which carries 
much lower complexity compared with the full SMS model. 

Proposition 2 (Euler-Lagrange Equations for Piecewise Constant SMS) The Euler-Lagrange 
equations for E[P,m \ I] in my\) are given by 



InP^ ' (25) 
-ISeApi + 2e-^pi{l - p^){l - 2pi) = {V) - X{mi - if, i = l:K, 

with Neumann conditions for all the ownerships pi 's along dfl. 



5.2 Computation of the Euler-Lagrange Equations 

Computationally, as well practiced in multivariate optimization problems, 1)241) and (|25() can be 
solved via the algorithm of alternating minimization (AM) |121 136j . The AM algorithm is closely 
connected to the celebrated EM (expectation-maximization) algorithm in statistical estimation 
problems with hidden variables |2()| I26j. In the current context, the ownership distributions pj's 
could be treated as the hidden variables. 

Like EM, the AM algorithm is progressive. Given the current {t = n) best estimation of the 
patterns U"" = {uf \ i = 1 : K), by solving 

P*" = argminp £^[P | LT'', /], 
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or equivalently, 

-18eApi + 2e-V(l-P*)(l-2pi) = (F")-A«-/)2, i = l:K, (26) 

with Neumann boundary conditions, one obtains the current best estimation of the ownerships 
P" = {p^ \ i = 1 : K). Subsequently, based on P", by solving 

LT'^+i = argmin^[J7 | P",/], 

or equivalently, 

-aAn, + {\p'^)ui = iXp^)I, i = 1 : K, (27) 

with Neumann boundary conditions, one completes a single round of pattern updating [/"" U"^^^. 
The same procedure applies to the piecewise constant soft Mumford-Shah equations in ((2SJ. 

Since the system ()27() is linear and decoupled, the main computational complexity resides in the 
integration of (|26|) . which is coupled and nonlinear due to the simplex constraint and the double- 
well potential in the energy. Define ej(x) = {ui{x) — /(x))^ and e = (ej | i = 1 : K). In order to 
solve 

-ISeApi + 2e-ipi(l - pi){l - 2pi) = {V) - Ae, (28) 
given e and y = F (P, U) = V{P, e) (see (EH), first notice that 

K 

K 



1 ^ 

(F(P,e)) = ^(-l8eAp, + ei + 2e-V(l-Pi)(l-2pi)) 



i=l 

i=l 1=1 

since Ylf=i Pi = ^ ^^^d A (^^^^ Pi^ ~ ^P^^^ double-potential force in (^5]) by 

Pi{l-Pi){l - 2pi) = pi{l - pif - pj{l - Pi). 
In combination, the nonlinear equation H28|) can then be solved iteratively: 

, p{j) _^ p{j+^) ^ . . . 

by the following linearization procedure: 



f!f) = -e, + (y(P<^-),e)) + 2e~' fpf)' fl - 



with Neumann adiabatic boundary conditions for all the channels i = 1 : K. This system of linear 
Poisson equations can be conveniently integrated using any elliptic solvers. The detailed numerical 
analysis on the convergence rates of the entire algorithm above, however, is still an open problem 
and well deserves some systematic investigation. 
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6 Computational Examples 



In this section, we present the computational results of the proposed soft Mumford-Shah model. 
Notice that the extension of the above SMS models to color images is straightforward by having the 
gray values iij's replaced by RGB vectors. (We, however, must remind the reader that perceptually 
RGB may not be the most ideal representation of colors compared with other nonlinear approaches, 
e.g., brightness-chromaticity [S] and HSV _5 .) 

Figures El and 0] illustrate the performance of the SMS model on two synthetic images with 
multiple phases. Figure 01 shows a typical T-junction and Figure a 3-phase image with a narrow 
bottleneck. Plotted in the figures are the hard segments obtained from the SMS model via the 
hardening formulae ((SJ and (jSJ. 



Figure too large for arXiv 



For the figure and the complete manuscript. 



go to: www.math.umn.edu/~jhshen 



Figure 3: Synthetic image of a T-junction: hard segmentation from the SMS model via "hardening" 
formulae ((JJ and (jSJ- The 120-degree regularization behavior at the junction point is also well known 
in the classical Mumford-Shah model ^33^. 

Plotted in Figure [51 are the hardened segments obtained from the soft Mumford-Shah segmen- 
tation model via formulae © and ©. For this application, a user specifies three small patches 
(three rectangles in this example) Qi,Q2, and Q3, and the SMS model proceeds with the extra 
interpolation conditions in (|14|) for the ownerships. 

In Figure O another example of a natural image is segmented via the SMS model and the 
"hardening" formulae ^ and (jSJ. A user supervises with three patches Qi,Q2 and Q3, and 
designates the two on the body to a pattern ownership Pbody and the third (from the ocean) to 
Pocean- If the three are treated as distinct patterns, the SMS model still works, but one needs an 
extra step of high-level vision processing (e.g., based on Grenander's graph models ^HI) to group 
the skin-tone and the purple-shirt patterns in order to capture the entire body faithfully. 

Finally, plotted in Figures ((Tj) and © are the ownerships from the SMS model based on the 
3-phase and 4-phase supervision separately in Figure El The stochastic nature of the outcomes 
(i.e., the softly transiting ownerships pj's instead of hard segmentation) is closer to the way how 
a human subject may perceive such a natural scene. In particular, the SMS model seems to be 
consistent with the most recent theory that hard pattern segments may not be absolutely necessary 
for natural scene recognition j24l I23j . 
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Figure too large for arXiv 



For the figure and the complete manuscript. 



go to: www.math.umn.edu/~jhshen 



Figure 4: Synthetic image of a narrow bottleneck: hard segmentation from the SMS model via 
"hardening" formulae © and (jSJ- The thickening regularization at the bottleneck junction can 
be explained similarly by the classical Mumford-Shah model for which minimum-surface or "soap- 
foam" behavior arises due to the surface tension energy. Also see the recent work by Kohn and 
Slastikov 22^ for the singularity analysis of a similar problem arising from micromagnetism. 
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